########### MANUSCRIPT: The Role of Education and Attitudes in Cooking Fuel Choice: Evidence from two states in India
########### JOURNAL: ENERGY FOR SUSTAINABLE DEVELOPMENT
########### AUTHORS: CARLOS F. GOULD/1 AND JOHANNES URPELAINEN/2
########### AFFILIATIONS: 1/COLUMBIA UNIVERSITY MAILMAN SCHOOL OF PUBLIC HEALTH AND 2/JOHNS HOPKINS SCHOOL OF ADVANCED INTERNATIONAL STUDIES
########### PURPOSE: THIS IS CODE #3 THAT MAKES DESCRIPTIVE FIGURES




# FIGURE 3 -------------------------------------------------------------

# Correlation matrix of Perception Indices and LPG ownership in (A) Kerala and (B)
# Rajasthan households using the Pearson's product-moment correlation.

# Note: there was some formatting in a photo editor of the lables and titles

### Kerala

k_perceptions_lpg_indices <- kerala[,c("Perception_Index_FirewoodValuable","Perception_Index_ChulhaDifficult",
                                     "Perception_Index_FirewoodUnavailableDifficult","Perception_Index_CookingHealth",
                                     "Perception_Index_GeneralHealth","Perception_Index_LPGGoodAffordable","Owns_LPG")]
colnames(k_perceptions_lpg_indices) <- c("Firewood is valuable", "Chulha is difficult to use", 
                                       "Firewood is unavailable, inconvenient",  "Cooking fuel type is related to health",
                                       "Generally progressive health norms", "LPG is a good cooking fuel and affordable", 
                                       "Has LPG")
k_perceptions_lpg_indices_matrix <- as.matrix(k_perceptions_lpg_indices)
k_perceptions_lpg_indices_corrs <- cor(k_perceptions_lpg_indices_matrix, use="na.or.complete", method = "pearson")  


k_corr <- GGally::ggcorr(data=NULL, cor_matrix=k_perceptions_lpg_indices_corrs, label=T, low="darkred", mid="white", high="steelblue",
       legend.size=16, size=6, label_size=8, label_round=2, layout.exp = 1, color="grey50", hjust=0.9)


### Rajasthan

r_perceptionslpg_indices <- rajasthan[,c("Perception_Index_LPGGood",
                                             "Perception_Index_LPGaffordable",
                                             "Perception_Index_CookingHealth",
                                             "Perception_Index_FirewoodValuable",
                                             "Perception_Index_FirewoodUnavailableDifficult",
                                             "Owns_LPG")]
colnames(r_perceptionslpg_indices) <- c("LPG is a good quality cooking and beneficial", 
                                        "LPG is desirable, but expensive", 
                                        "Cooking fuel type is related to health", 
                                        "Chulha and firewood are valuable",
                                        "Firewood is unavailable and inconvenient", 
                                        "Has LPG")
r_perceptionslpg_indices_matrix <- as.matrix(r_perceptionslpg_indices)
r_perceptionslpg_indices_corrs <- cor(r_perceptionslpg_indices_matrix, use="na.or.complete")  

r_corr <- ggcorr(data=NULL, cor_matrix=r_perceptionslpg_indices_corrs, label=T, low="darkred", mid="white", high="steelblue",
       legend.size=16, size=6, label_size=8, label_round=2, layout.exp = 1, color="grey50", hjust=0.9)


cowplot::plot_grid(
  
  k_corr, r_corr, 
  labels = c("(A) Kerala Perceptions Indices and LPG", "(B) Rajasthan Perception Indices and LPG"),
  nrow = 2
  
)


# FIGURE A1 -------------------------------------------------------------
# Directed acyclic graph (DAG) outlining a priori causal theory of the relationship of
# education with LPG adoption.
# Made by hand in powerpoint



# FIGURE A2 -------------------------------------------------------------
# Correlation matrix between control variables demonstrating little association among
# covariates.

corr_vars <- merged_regs[,c("Exp_log", 
                            "NumberAdults", "NumberChildrenUnder18",
                            "AgeRespondent",
                            "Caste_NotGeneral",
                            "Religion_NonHindu",
                            "Rural")]


corr_vars_matrix <- as.matrix(corr_vars)
corr_vars_matrix_corrs <- cor(corr_vars_matrix, use="na.or.complete")  

corr_vars_plot <- ggcorr(data=NULL, cor_matrix=corr_vars_matrix_corrs, label=T, low="darkred", mid="white", high="steelblue",
                 legend.size=16, size=6, label_size=8, label_round=2, layout.exp = 1, color="grey50", hjust=0.9)

corr_vars_plot


# FIGURE A4 -------------------------------------------------------------
# Density plot showing the distribution of self-reported days to refill an LPG cylinder in
# Kerala and Rajasthan LPG-owning households.


ggplot(merged_regs, aes(x=q126_1)) + 
  geom_density(aes(y = ..scaled..), fill="black") + 
  theme_bw() + 
  xlab("Days before LPG cylinder refill") + 
  ylab("Density (scaled to 1)") + 
  theme(plot.title=element_text(size=24),
        plot.subtitle = element_text(size=24),
        axis.text = element_text(size=18),
        axis.title = element_text(size=22),
        legend.title = element_blank(),
        legend.position = c(0.85,0.95),
        strip.text = element_text(size=20)) +
  facet_grid(state~.)


# FIGURE A5 -------------------------------------------------------------
# Density plots showing the distribution of perception indices in LPG and non-LPG
# owning households in the Kerala and Rajasthan study samples.


# Kerala

kerala$Owns_LPG_labeled <- ifelse(kerala$Owns_LPG==1, "LPG", "No LPG")

p1_dist <- ggplot(kerala, aes(x=Perception_Index_FirewoodValuable_Norm)) + 
  geom_density(aes(y = ..scaled..), fill="black") + theme_bw() + 
  ggtitle("Chulha and woodfuel are valuable") +
  xlab("") + ylab("Density (scaled to 1)") + 
  theme(plot.title=element_text(size=24),
        plot.subtitle = element_text(size=24),
        axis.text = element_text(size=18),
        axis.title = element_text(size=22),
        legend.title = element_blank(),
        legend.position = c(0.85,0.95),
        legend.text = element_text(size=14)) + 
  facet_grid(Owns_LPG_labeled~.)

p2_dist <- ggplot(kerala, aes(x=Perception_Index_ChulhaDifficult_Norm)) + 
  geom_density(aes(y = ..scaled..), fill="black") + theme_bw() + 
  ggtitle("Using the chulha is difficult") +
  xlab("") + ylab("") + 
  theme(plot.title=element_text(size=24),
        plot.subtitle = element_text(size=24),
        axis.text = element_text(size=18),
        axis.title = element_text(size=22),
        legend.title = element_blank(),
        legend.position = c(0.85,0.95),
        legend.text = element_text(size=14)) + 
  facet_grid(Owns_LPG_labeled~.)

p3_dist <- ggplot(kerala, aes(x=Perception_Index_FirewoodUnavailableDifficult_Norm)) + 
  geom_density(aes(y = ..scaled..), fill="black") + theme_bw() + 
  ggtitle("Woodfuels are not available, inconvenient to use") +
  xlab("") + ylab("Density (scaled to 1)") + 
  theme(plot.title=element_text(size=20),
        plot.subtitle = element_text(size=24),
        axis.text = element_text(size=18),
        axis.title = element_text(size=22),
        legend.title = element_blank(),
        legend.position = c(0.85,0.95),
        legend.text = element_text(size=14)) + 
  facet_grid(Owns_LPG_labeled~.)

p4_dist <- ggplot(kerala, aes(x=Perception_Index_CookingHealth_Norm)) + 
  geom_density(aes(y = ..scaled..), fill="black") + theme_bw() + 
  ggtitle("Cooking fuel choices affect health") +
  xlab("") + ylab("") + 
  theme(plot.title=element_text(size=24),
        plot.subtitle = element_text(size=24),
        axis.text = element_text(size=18),
        axis.title = element_text(size=22),
        legend.title = element_blank(),
        legend.position = c(0.85,0.95),
        legend.text = element_text(size=14)) + 
  facet_grid(Owns_LPG_labeled~.)

p5_dist <- ggplot(kerala, aes(x=Perception_Index_GeneralHealth_Norm)) + 
  geom_density(aes(y = ..scaled..), fill="black") + theme_bw() + 
  ggtitle("General health attitudes are progressive") +
  xlab("") + ylab("Density (scaled to 1)") + 
  theme(plot.title=element_text(size=24),
        plot.subtitle = element_text(size=24),
        axis.text = element_text(size=18),
        axis.title = element_text(size=22),
        legend.title = element_blank(),
        legend.position = c(0.85,0.95),
        legend.text = element_text(size=14)) + 
  facet_grid(Owns_LPG_labeled~.)

p6_dist <- ggplot(kerala, aes(x=Perception_Index_LPGGoodAffordable_Norm)) + 
  geom_density(aes(y = ..scaled..), fill="black") + theme_bw() + 
  ggtitle("LPG is good and affordable") +
  xlab("") + ylab("") + 
  theme(plot.title=element_text(size=24),
        plot.subtitle = element_text(size=24),
        axis.text = element_text(size=18),
        axis.title = element_text(size=22),
        legend.title = element_blank(),
        legend.position = c(0.85,0.95),
        legend.text = element_text(size=14)) + 
  facet_grid(Owns_LPG_labeled~.)

grid.arrange(p1_dist, p2_dist, p3_dist, p6_dist, p5_dist, p4_dist)

k_p_dist <- grid.arrange(p1_dist, p2_dist, p3_dist, p6_dist, p5_dist, p4_dist)



# RAJASTHAN

rajasthan$Owns_LPG_labeled <- ifelse(rajasthan$Owns_LPG==1, "LPG", "No LPG")

p1_dist <- ggplot(rajasthan, aes(x=Perception_Index_FirewoodValuable)) + 
  geom_density(aes(y = ..scaled..), fill="black") + theme_bw() + 
  ggtitle("Chulha and woodfuel are valuable") +
  xlab("") + ylab("Density (scaled to 1)") + 
  theme(plot.title=element_text(size=24),
        plot.subtitle = element_text(size=24),
        axis.text = element_text(size=18),
        axis.title = element_text(size=22),
        legend.title = element_blank(),
        legend.position = c(0.85,0.95),
        legend.text = element_text(size=14)) + 
  facet_grid(Owns_LPG_labeled~.)

p2_dist <- ggplot(rajasthan, aes(x=Perception_Index_FirewoodUnavailableDifficult)) + 
  geom_density(aes(y = ..scaled..), fill="black") + theme_bw() + 
  ggtitle("Woodfuels are not available, inconvenient to use") +
  xlab("") + ylab("") + 
  theme(plot.title=element_text(size=20),
        plot.subtitle = element_text(size=24),
        axis.text = element_text(size=18),
        axis.title = element_text(size=22),
        legend.title = element_blank(),
        legend.position = c(0.85,0.95),
        legend.text = element_text(size=14)) + 
  facet_grid(Owns_LPG_labeled~.)

p3_dist <- ggplot(rajasthan, aes(x=Perception_Index_CookingHealth)) + 
  geom_density(aes(y = ..scaled..), fill="black") + theme_bw() + 
  ggtitle("Cooking fuel choices affect my health") +
  xlab("") + ylab("Density (scaled to 1)") + 
  theme(plot.title=element_text(size=20),
        plot.subtitle = element_text(size=24),
        axis.text = element_text(size=18),
        axis.title = element_text(size=22),
        legend.title = element_blank(),
        legend.position = c(0.85,0.95),
        legend.text = element_text(size=14)) + 
  facet_grid(Owns_LPG_labeled~.)

p4_dist <- ggplot(rajasthan, aes(x=Perception_Index_LPGaffordable)) + 
  geom_density(aes(y = ..scaled..), fill="black") + theme_bw() + 
  ggtitle("LPG is desirable, but expensive") +
  xlab("") + ylab("") + 
  theme(plot.title=element_text(size=24),
        plot.subtitle = element_text(size=24),
        axis.text = element_text(size=18),
        axis.title = element_text(size=22),
        legend.title = element_blank(),
        legend.position = c(0.85,0.95),
        legend.text = element_text(size=14)) + 
  facet_grid(Owns_LPG_labeled~.)

p5_dist <- ggplot(rajasthan, aes(x=Perception_Index_LPGGood)) + 
  geom_density(aes(y = ..scaled..), fill="black") + theme_bw() + 
  ggtitle("LPG is good quality cooking and beneficial") +
  xlab("") + ylab("Density (scaled to 1)") + 
  theme(plot.title=element_text(size=24),
        plot.subtitle = element_text(size=24),
        axis.text = element_text(size=18),
        axis.title = element_text(size=22),
        legend.title = element_blank(),
        legend.position = c(0.85,0.95),
        legend.text = element_text(size=14)) + 
  facet_grid(Owns_LPG_labeled~.)


grid.arrange(p1_dist, p2_dist, p3_dist, p4_dist, p5_dist)

r_p_dist <- grid.arrange(p1_dist, p2_dist, p3_dist, p4_dist, p5_dist)

plot_grid(
  
  k_p_dist, r_p_dist,
  labels = c("(A) Kerala Perceptions Indices", "(B) Rajasthan Perceptions Indices"),
  nrow = 2
  
)

# FIGURE A6 -------------------------------------------------------------
# Density plots showing the distribution of the price (in rupees) paid by households for
# their most recent LPG cylinder purchase in the Kerala and Rajasthan study samples. Prices above
# 1,500 rupees were removed due to implausibility.

ggplot(merged_regs, aes(x=q127a_1)) + 
  geom_density(aes(y = ..scaled..), fill = "black") + 
  theme_bw() + 
  xlim(c(0, 1500)) + 
  xlab("Last price paid for LPG cylinder") + 
  ylab("Density (scaled to 1)") + 
  theme(plot.title=element_text(size=24),
        plot.subtitle = element_text(size=24),
        axis.text = element_text(size=18),
        axis.title = element_text(size=22),
        legend.title = element_blank(),
        legend.position = c(0.85,0.95),
        strip.text = element_text(size=20)) +
  facet_grid(state~.)
